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This is Part II of the paper series on data-compatible T-matrix completion (DCTMC), which is 
a method for solving nonlinear inverse problems. Part I of the series contains theory and here we 
present simulations for inverse scattering of scalar waves. The underlying mathematical model is the 
scalar wave equation and the object function that is reconstructed is the medium susceptibility. The 
simulations are relevant to ultrasound tomographic imaging and seismic tomography. It is shown 
that DCTMC is a viable method for solving strongly nonlinear inverse problems with large data 
sets. It provides not only the overall shape of the object but the quantitative contrast, which can 
correspond, for instance, to the variable speed of sound in the imaged medium. 


I. INTRODUCTION 

This paper is Part II of the series on solving nonlinear 
inverse scattering problems (ISPs) by data-compatible 
T-matrix completion (DCTMC). Part I [l| contains the¬ 
ory and here we report initial numerical simulations for 
the three-dimensional diffraction tomography with scalar 
waves. This problem arises, in particular, in tomographic 
ultrasound imaging ii and in seismology 0-0 ■ We 
note that DCTMC can be applied in a very similar man¬ 
ner to the problems of inverse electroma gne tic scatter¬ 
ing and diffuse optical tomography |l^ . [l^ . How¬ 

ever, in this paper we probe the medium with scalar prop¬ 
agating waves. By doing so, we avoid, on one hand, the 
additional complexity related to the vectorial nature of 
electromagnetic fields and, on the other hand, the severe 
ill-posedness of the ISP associated with the exponential 
decay of diffuse waves. Therefore, we can focus on the 
effects of nonlinearity of the ISP. For the same reason, 
no noise will be added to the simulated data. 

We have studied convergence and computational per¬ 
formance of several modifications of the DCTMC algo¬ 
rithm. We start with the streamlined iteration cycle with 
both computational shortcuts as described in [l| used. 
This approach minimizes computational time per one it¬ 
eration of DCTMC. However, it is not the most efficient 
algorithm in terms of the number of iterations required 
for convergence. We then implement several modifica¬ 
tions and improvements of DCTMC. In particular, we 
utilize weighted summation to the diagonal for the “force- 
diagonalization operator” D[-], which is used to reduce 
the iteratively-updated interaction matrix 14 to a diag¬ 
onal matrix D^. This makes the use of Computational 
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Shortcut 2 impossible. However, weighted summation to 
the diagonal uses explicitly the off-diagonal elements of 
14 and is, therefore, more in line with the main idea of 
DCTMC. Correspondingly, the convergence of DCTMC 
is significantly improved. With all the improvements, the 
number of required iterations in a representative test case 
was reduced from 900 (for the standard algorithm) to 75. 
We also demonstrate how physical constraints and checks 
for sparsity can be embedded in the DCTMC algorithm. 
Application of these checks can be viewed as a method 
of regularizing the nonlinear ISP. 

The remainder of this paper is organized as follows. 
In Sec. im we describe the procedure for discretization 
of the scalar wave equation. Technical details of the nu¬ 
merical procedures, description of the targets and of the 
source-detector arrangements used are given in Sec. IHII 
Numerical results for the DCTMC algorithm without any 
modifications are shown in Sec. El In Sec. El we illus¬ 
trate several methods to improve the rate of convergence 
of DCTMC. In Sec. IVII we show that DCTMC can pro¬ 
vide good results when the Gauss-Newton method fails. 
Finally, Sec. IVHI contains a brief discussion of the ob¬ 
tained results. 


II. DISCRETIZATION OF THE SCALAR WAVE 
EQUATION 

Consider a scalar field u(r) (e.g., the pressure wave in 
ultrasound imaging) and the wave equation 

[V^ -I- A:^e(r)] M(r) = —47rfc^(7(r) , (1) 

where q{T) is the source function, e(r) = 1 outside of the 
bounded domain D (the sample) and the factor —47rfc^ in 
the right-hand side has been introduced for convenience. 
We work in the frequency domain and assume that the 
wave number k = uj/c is fixed, where c is the velocity of 
waves in free space, that is, outside of D. 
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Since we wish to implement a numerical procedure 
for reconstructing e(r) from external measurements, the 
problem must be suitably discretized and the unknown 
function must be represented by a finite number of de¬ 
grees of freedom (voxels). To this end, we follow the gen¬ 
eral approach of the discrete dipole approximation that 
was originally developed for Maxwell’s equations 
A related method for the scalar diffusion equation was 
described by us in [l^ . 

We start by re-writing m identically as 

(V^-I-u(r) =-47rfc^ [x(r)u(r)-I-q(r)] , (2) 


where 


x(i-) 


e(i~) - 1 
47r 


( 3 ) 


This quantity can be referred to as the susceptibility 
of the medium. Inverting the differential operator in 
the left-hand side of ([2]), we obtain Lippmann-Schwinger 
equation 

u(r) = Minc(r)-I- J Go{r,r')x{r')u{r')(fr', (4) 

where 


= (!■) 


= J Go{r, r')q{r')(fir' 


( 5 ) 


is the incident field and 

Go(r,r') = k 


2 exp {ik\r — r'|) 


( 6 ) 


is the free-space Green’s function of the wave equation, 
which satisfies 


where 


Cn — r/,inc(l’n) ■ (H) 

Equation (fTUl) can be reasonably expected to have a so¬ 
lution, which then gives a discrete approximation to the 
continuous field it(r). 

In principle, we can compute the integrals in the right- 
hand side of OT analytically or numerically with any 
degree of precision. However, doing so is not practically 
useful because there is already an approximation involved 
in writing m- Therefore, the conventional approach is 
to approximate the integrals by expressions that are easy 
to compute if Gq is known analytically. To obtain such 
expressions, we consider the two cases n ^ m and n = m 
separately. 

For n ^ m, the standard approximation is 



Go(r„,r)d^r 


h Go(i‘^,r,72) , 


n ^ m . (12) 


This approximation is not very good for neighboring cells 
but generally believed to be adequate for cubic lattices (it 
would be precise for spherical region). However, it should 
not be applied unaltered for orthorhombic lattices with 
primitive lattices of non-equal length [T^. 

For n = m, the integrand contains a singularity and 
a more careful evaluation of the integral is required. It 
is standard to assume that fc/i <C 1 (otherwise, the dis¬ 
cretization is not a valid approximation to the underlying 
differential equation) so that the exponent in (jS]) can be 
approximated as 


exp (iA:|r — r'l) Ri 1-I-jfc|r — r'l . (13) 


(V^ -I- k'^) Go(r, r') = —47rfc^5(r — r') . (7) 

Note that the second term in the right-hand side of (|4]) 
is the scattered field Mscatt(i’), that is, 

Mscatt(r) = j Go{r,r')x{r')u{r')(fr' . (8) 

The total field is given by a sum of the incident and 
scattered components, u = Ui„c + Ugcatt- 

We now proceed with discretization of the integral 
equation Let the sample be rectangular and dis¬ 

cretized into cubic voxels C„ (n = 1 ,..., IVy) of volume 
each, where Ny is the total number of voxels. We then 
make the following approximation [II: 

X(r) = Xn and u(r) = IF r S C„ . (9) 

Here Xn and are constants. Obviously, with this ap¬ 
proximation used, o can not hold for all values of r. 
However, by restricting attention to the points r = r„, 
where r„ are the centers of the respective voxels C„, we 
obtain a set of algebraic equations 

Ny 

4" ^ ^ Xm^m / Go(l‘m I’)(i V , (1^) 


We then obtain 



Go(r„,r)d^r 



X 




x'^ + y'^ + z'^ 



{khfy + ikh) , (14) 


where 

C = log(26-f 15\/3) -7r/2a!2.38 . (15) 

The imaginary part of the right-hand side of (fCT) is the 
first non-vanishing radiative correction to the voxel po¬ 
larizability. As is known for the electromagnetic case, ac¬ 
counting for this correction is important to ensure energy 
conservation of the scattering process Note that 

making higher-order approximations in (1131) is known to 
produce higher-order dynamic corrections. For the elec¬ 
tromagnetic case, the relevant results are given in [2ll . [^ . 
Also, an alternative estimate of the parameter ^ can be 
obtained if we replace integration over a cubic voxel by 
integration over a ball of equal volume. This approach 
results in a much simpler integral and does not affect the 
radiative correction; however, it gives a slightly different 
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value for namely, ^ = (97r/2)^/^ Ri 2.42. The two ap¬ 
proaches to computing ^ can be applicable to different 
physical situations. For example, integration over a ball 
may be more appropriate if we wish to describe scat¬ 
tering by a collection of spherical particles rather than 
voxelization of a medium on a cubic grid. A review of 
different definitions of the parameter ^ for the electro¬ 
magnetic case and a summary of several relevant expres¬ 
sions, including (ITSl) . are given in (^ . 

Further, it is convenient to introduce the ’’moments” 

dn = h^XnUn ■ (16) 

In terms of dn, and with the account of the integration 
results dm and (HI, the system of equations (Uni takes 
the form 




where 


dn — O^n I Cyi -f ^ ( r nmdn 

\ m—1 

h^Xn 


— 


1 - -I- ikh)xn 

is the polarizability of n-th voxel and 

^nm (f ^nm) ^0 ^m) 


(17) 


(18) 


(19) 


is the interaction matrix. Note that Fnm has zeros on 
the diagonal. 

If all polarizabilities are known, then equation 
(Hzl) is the statement of the forward scattering problem. 
Namely, given the incident field in all voxels, e„, find the 
moments dn by solving (II 3 ; then use the formula 

Nv 

'ascatt(^^d) — 'y G(}(^V(i,Tn)dn (^6) 

n—1 

[a discretized version of (|H])] to find the scattered field 
at an arbitrary point of observation r^;. Obviously, the 
forward problem is linear with respect to the unknowns, 
dn • 

The inverse problem is stated differently and is, in gen¬ 
eral, nonlinear. Our goal is to use multiple measurements 
of the scattered field Mscatt(rd) due to multiple external 
point sources of the form g(r) = (5(r — Ts), where and 
r^, are the positions of the detector and the source, to 
find all the voxel polarizabilities a„. In the matrix form, 
the ISP is stated as follows. Let us define the interaction 
matrix V as the Ny x Ny matrix that contains on the 
diagonal and zeros elsewhere. Then we can re-write m 
in matrix notations as 


\d) = V{\e) + r\d)) , (21) 

where jd) and |e) are the vectors of the length Ny with 
the elements dn and e„, respectively. The formal solution 
to (1^ is 


( 22 ) 


where 


T[V] = {I-Vr)-^V (23) 


is the T-matrix. Let us assume that we measure the 
scattered field at a set of points Vdk £ 'Fd, k = 1,... ,Nd, 
and let fk = itscatt(i’dfc)- Then |/) is a vector of the 
length Nd and we can use o to write |/) = A\d), where 
Akn = Go(i'dfe,i'n)- Similarly, we have |e) = B\q), where 
Bni = Go{rn, r^i) and r*/ e E*, ^ = 1 ,..., is the set of 
points from which we illuminate the medium. By turning 
on one source at a time and by measuring the scattered 
field at all points Tdk for each source, we measure all 
elements of the data matrix (l>ki. As follows from (I22L 
the data matrix satishes the equation 


AT[V]B = $ . (24) 


This is equivalent to Eqs. 4 and 8 of [T|. As in [l|, the 
matrices A and B are obtained directly by sampling the 
Green’s function Go(r,r'), for which we have given a 
specific expression dS) that corresponds to the physical 
problem considered. The only fine (and slightly nontriv¬ 
ial) point here is the way in which the matrix F was 
obtained. For n to, we still have Fnm = Go{rn,Vm), 
i.e., the off-diagonal matrix elements were obtained by 
sampling Gq. However, the diagonal elements of F can 
not be obtained by straightforward sampling due to the 
singularity of Gq. The discretization procedure described 
above results in Fnn = 0 provided that we use the renor¬ 
malized polarizabilities (1181) to quantify the response 
of individual voxels to external excitation. 

The above approach to obtaining the renormalized po¬ 
larizabilities is rather standard in physics and goes back 
to the Lorentz’s idea of local-field correction (in electro¬ 
magnetics). However, it is useful to keep in mind that 
the local-held correction and the related renormalization 
are mathematically related to the singularity of the free- 
space Green’s function. From the purely algebraic point 
of view, we can point out that the renormalization in 
question is a special case of the identical transformation 


(J - FT)-V = {I- v'F'yW 


where V = PV , F' = F — V~^{I — P~^) and P is any 
invertible matrix (although V~^ appears in the above 
expression, invertibility of V is not really required). In 
particular, if F has the diagonal part D, then we can take 
P = {I — VD)~^ (assuming that this inverse exists) and 
the renormalized matrix F' will have a zero diagonal. It 
also should be pointed out that the renormalized polariz¬ 
abilities derived here are adequate for cubic lattices but 
require modihcation in the case of orthorhombic lattices 
with primitive vectors of unequal length [l^ . 


\d) = {I-VF)-^V\e)=T[V]\e) 



III. DETAILS OF THE NUMERICAL 
ALGORITHMS USED 


B. Application of physical constraints 
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A. Definition of the unknowns 

The above inverse problem has been formulated in 
terms of the unknown voxel polarizabilities However, 
the physical quantity of interest is the voxel permittivity 
e„ or the susceptibility Xn- To relate the latter to the 
former, we can use m and (|T8|) to obtain in a straight¬ 
forward manner 


__ 

1 -I- ikh){an/h?) 


(25) 


The permittivities can be trivially related to Xn by 
using and we will focus on the latter quantities in 
the remainder of this paper. 

Equation (EH) is an analog of the Maxwell Garnett 
formula written for scalar waves (with the account of the 
first non-vanishing radiative correction) and its inverse is 
the Clausius-Mossotti relation. It can be seen that the 
relation between and a„ is itself nonlinear. In other 
words, we have removed some nonlinearity of the ISP an¬ 
alytically. The nonlinearity in question is solely due to 
the self-interaction of a given voxel and it was accounted 
for by the procedure of renormalization of the polariz¬ 
ability. In other words, we can view the numerator h^Xn 
in the right-hand side of (flSl) as the bare polarizability 
and the complete expression as the renormalized polariz¬ 
ability. However, the nonlinearity of the ISP that is due 
to the interaction of different voxels can not be removed 
so simply and we will have to deal with it numerically. 

Still, the above discussion gives some validity (al¬ 
though not a rigorous proof) to the proposition that it 
is advantageous to formulate the ISP in terms of the po¬ 
larizabilities an rather than in terms of the susceptibil¬ 
ities Xn- Indeed, consider the case when there is only 
one voxel and we wish to determine its properties (either 
a or x) from external measurements by means of some 
generic iterative scheme that is applicable to a more gen¬ 
eral problem (e.g., involving many interacting voxels). If 
we view a as the fundamental unknown, we arrive at a 
well-posed linear equation of the form Aa = b {A ^ 0 
and b are numbers), which can be solved in just one it¬ 
eration. However, if we view y as the fundamental un¬ 
known, then we will be solving iteratively an equation 
of the form ^x/(l — Px) = b, which can take several 
iterations depending on the value of the coefficient /3. 

In the simulations reported below, we formulate the 
model and display the reconstructions in terms of the 
susceptibilities Xn, which we assume to be real-valued. 
However, to generate the data function, we compute the 
set of a„’s according to (fT51) . Then we ’’pretend” that 
a^s are unknown and, viewing these quantities as the 
fundamental unknowns, solve the ISP. We then convert 
the reconstructed values of a„’s to XnS by using (1^ and 
display the latter quantities in all figures. 


As discussed in [l[ , the iterative procedure of DCTMC 
can benefit from applying physical constraints to the 
unknowns as a form of regularization. The constraints 
are usually based on general physical principles but can 
also accommodate some a priori knowledge about the 
medium. In this work, only fundamental constraints are 
used. Specifically, if the medium is passive (that is, non¬ 
amplifying), then we know that Imxn >0. If it is also 
known that the medium is transparent (non-absorbing), 
then Imx„ = 0. If we view the polarizabilities a„ as the 
fundamental unknowns in the iterative DCTMC process, 
we can apply this physical constraint in the following 
manner. We first notice that 



'Imxn 

JXnP 


+ {khf 


< -{khf . 


(26) 


Let’s say, a numerical iteration of DCTMC has produced 
a set of a„. To enforce the condition (l26t . we can apply 
the following transformation: 


Re(l/a„) 


1 

i max [—Im(l/Q:„), k^] 


(27) 


If, in addition, we know that the sample is non-absorbing 
so that Imx„ = 0 for all voxels and strict equality in the 
last relation in (1261) holds, then we can apply an even 
sharper constraint by writing 


1 

Re(l/a„) — ik^ 


(28) 


In the simulations reported below, we have assumed that 
the sample is transparent and used the constraint (l28l) at 
each iteration of DCTMC. We note that in many cases 
considered, DCTMC produces very similar results even 
without applying the constraint, but then the conver¬ 
gence is somewhat slower. The only instances in which 
we found that the physical constraint is critical were the 
cases with very strong nonlinearity. 


C. Account of sparsity 

In many practical cases we can expect the target to be 
in some sense sparse. This means that many of the voxels 
have zero (or relatively small) susceptibilities Xn, but we 
do not know a priori where these voxels are located or 
how many such voxels exist in the computational domain. 
We will refer to such voxels as ’’noninteracting” as the 
moments ci„ of these voxels are relatively small, do not 
interact effectively with the moments of other voxels, and 
produce a negligible input to the scattered field. 

In DCTMC, it is possible to take the sparsity into 
account in an adaptive manner without actually know¬ 
ing whether the target is sparse or not or where the 
non-interacting voxels are located. In the simulations 
of Sec. HV] below, we have used the following rather ad 
hoc algorithm: 
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1: Run 50 iterations normally. 

2: Then every 20 iterations check whether some sus¬ 
ceptibilities Xn satisfy \xn\ < Xmax/100, where 
Xmax = max„ Ixnl- 

3: If a given voxel satisfies the above condition 3 
checks in a row, the corresponding Xn is set to zero. 

4: The voxels with zero Xn (as determined in the 
previous step) are declared to be non-interacting 
and are excluded from the computational domain. 
When this happens, we repeat the initial DCTMC 
setup, but now for a smaller number of interacting 
voxels Ny. This results in a smaller computational 
time per a subsequent iteration. 

5: The process is repeated with the following modifi¬ 
cations. After 200 iterations, checks are made every 
10 iterations, and after 400 iterations, the relative 
threshold for determining a non-interacting voxel is 
reduced to the factor of 60, and after 600 iterations 
the relative threshold is reduced to the factor of 40. 

Note that all integer constants used in the above al¬ 
gorithm are adjustable. In Sec. lYl we have used some¬ 
what different values of these constants (as noted in each 
particular case) and in some simulations we did not use 
sparsity checks at all. 

The procedure of selecting the noninteracting voxels 
can be described as iterative ’’roughening” of the tar¬ 
get. Any iterative numerical reconstruction is expected 
to produce small but nonzero values in the regions of the 
computational domain where the reconstructed function 
is really zero. Keeping these small values in subsequent 
iterations is nothing but a computational burden. The 
procedure described here removes this burden without 
affecting the final images in a significant way as long as 
the noninteracting voxels are not assigned incorrectly. In 
the simulations shown below, such incorrect assignment 
occurred in the strong nonlinearity regime. We note that 
such occurrences can be easily ’’diagnosed” by monitor¬ 
ing the error of the matrix equation (1241) as is explained 
below. If an incorrect assignment does occur, one can po¬ 
tentially alleviate the problem by adjusting the constants 
used by the algorithm or by not using sparsity checks at 
all. 

Finally note that the homogeneous background against 
which the roughening is performed can be different from 
zero. Thus, in optical tomography, it is conventional to 
assume that the background properties of the medium 
are known yet different from those of free space In 
this case, a more complicated Green’s function Gq must 
be used. We emphasize that Go can be computed ana¬ 
lytically for many regular geometries of the background 
medium. In this paper, we have taken the homogeneous 
background to be free space (vacuum) with the sole pur¬ 
pose of being able to use the mathematically-simple func¬ 
tion given in (jH]). 


D. Model targets and arrangement of sources and 
detectors 

We have used three kinds of targets in the simulations, 
which we refer to as large, small and tiny. A target of 
a particular kind has always the same ’’shape” but can 
have varying degrees of contrast. The contrast is defined 
in terms of the susceptibility x(i’)- Mathematically, this 
means that x(i’) for ^ given target can be written as 

X(r) = xo0(r) , (29) 

where 0 < 0(r) < 1 is the shape function (always the 
same for a target of a given kind) and xo > 0 is the vari¬ 
able amplitude. Obviously, the larger is xoj the stronger 
is the nonlinearity of the ISP. Note that all reconstruc¬ 
tions shown below display the ratio Xn/xoj which, ideally, 
should coincide with the shape function. However, we did 
not use any a priori knowledge about xo to obtain the 
reconstructions. The normalization to xo was performed 
a posteriori, which allowed us to use the same color scale 
in all figures. 

The majority of simulations shown below are per¬ 
formed for the small target, which is discretized on a 
16 X 16 X 9 grid {Ny = 2,304). DCTMC reconstruc¬ 
tions were also demonstrated for the large target, which 
is discretized on a 30 x 30 x 15 grid {Ny = 13, 500). The 
tiny target (8 x 8 x 4, = 256) is used only in Sec. lyTl 

below for the purpose of comparison of DCTMC with 
the Gauss-Newton methods since the latter did not yield 
useful results for the small target and were too compu¬ 
tationally inefficient to be applied for the large target. 

The shapes of all targets are shown in Fig. [TJ It can 
be seen that the small target contains two rectangular 
inclusions in a homogeneous zero background. One in¬ 
clusion is of the size 6 x 6 x 3 (in units of h) and has 
0 = 1.0 and the second inclusion is of the size 5x5x2 
and has 0 = 0.857. Note that the smaller and the larger 
inclusions touch at one corner. The large target contains 
one rectangular inclusion of the size 6x6x4 and with 
0 = 1 and another inclusion of the size 10 x 10 x 6 and 
with 0 = 0.75. These two inclusions do not quite touch. 
The tiny target contains an inhomogeneity in the form of 
two squares of the size 4 x 4 x 1 in located in the two in¬ 
terstitial slices, as is shown in the figure. The two square 
regions overlap in the central 2x2 area when viewed from 
top. In terms of the shape function, we have 0(r) = 1 
inside the inhomogeneity and 0(r) = 0.5 in the back¬ 
ground. Thus, all voxels in the tiny target are different 
from free space. 

It can be remarked that the tiny target does not rep¬ 
resent an accurate discretization of the underlying differ¬ 
ential equation. However, it poses a nonlinear ISP in its 
own right, and various inversion methods can be tested 
using this model. 

The relation between the discretization step h and the 
free-space wave number k is kh = 0.2 in all cases. We can 
use this numerical value to estimate roughly the degree 
of nonlinearity of the ISP. We generally expect the ISP 
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FIG. 1. The shapes of all three targets and the color scheme 
used in all figures below. Even though all model shape func¬ 
tions satisfy (by definition) the condition 0 < 0 < 1, the 
color scale can be used to represent any quantity in the 
range [—0.3,1.3]. The allowances are made because the re¬ 
constructed values of Xn/Xo can become negative or larger 
than unity. Reconstructed values that are smaller than —0.3 
will be shown by the uniform black color and values that are 
larger than 1.3 will be shown as white. 


to be very nonlinear when the phase shift between two 
waves - one propagating through an inhomogeneity (say, 
in the z-direction) and another propagating through the 
background medium - becomes of the order of 7r/2 or 
larger. This phase shift is given by 

= {kh)n (^^yi + 47rxo0* - l) , (30) 

where n is the inhomogeneity depth in units of h, xo is the 
contrast, and 0i is the value of the shape function inside 
the inhomogeneity. For example, the largest contrast we 
used for the small target is xo = 1-75. The shape func¬ 
tion inside the inhomogeneity that is n = 3 voxels deep 
is 0i = 0.857. This corresponds to Aip O.GGtt. The 
nonlinearity in this case is expected to be very strong. 
Of course, this analysis does not take into account multi¬ 
ple scattering between different inhomogeneities. In fact, 
we’ll see that the linearized reconstructions break down 
for a much smaller contrast, i.e., at xo ~ 0.1. 

We now turn to the source-detector arrangements. 

In the case of the small target, the mesh of sources is 
a 22 X 22 rectangular grid with the same spacing h as 
was used to discretize the sample. The plane of sources 
is centered symmetrically near one of the 16 x 16 faces of 
the sample so that there are three rows and columns of 


sources extending past the sample surface in each direc¬ 
tion. The plane of sources is displaced by h/2 from the 
physical surface of the sample. The mesh of detectors is 
identical in dimension and located on the other side of 
the sample. Regarding the size of the data set, we have 
Ns = Nd = 484 and the total number of data points 
is NgNd = 234, 256. Note that the additional rows and 
columns of sources and detectors extending past the sam¬ 
ple surface are needed to achieve the best possible linear 
reconstruction in the weak nonlinearity regime. Adding 
even more rows and columns into the source/detector 
meshes does not improve this result any further. 

In the case of the large target, the grids of sources and 
detectors are of the size 38 x 38 so that four rows and 
columns of sources or detectors are extending past the 
sample surface in each direction. The source/detector 
meshes are displaced from each of the two 30 x 30 faces 
of the large target by h/2, as is also the case for the 
small target. The data set size is defined by the following 
numbers: Ns = Nd = 1,444 and NsNd = 2,085,136. 

For the tiny target, the grids of sources and detectors 
were of the size 10 x 10, {Nd = Ng = 100) and centered 
about the sample on all sides. 


E. Methods for linearized reconstruction 


Applying a nonlinear solver to an ISP is meaningful 
only if the linearized solution to the same problem has 
failed. Therefore, we compare the results of nonlinear 
reconstruction to those obtained by linearized inversion. 
We have used two different approaches to this problem. 

First, we have used the traditional approach of com¬ 
puting the pseudoinverse of the matrix K (see Remarks 
2 and 3 in i)- This approach leads to a linear prob¬ 
lem whose size grows very rapidly with the size of the 
data set and the number of voxels. While still feasi¬ 
ble in the case of the small target, application of this 
method for the large target is already problematic. In¬ 
deed, the matrix K for the large target has the total 
number of elements NgNdNy = 28,149,336,000. Stor¬ 
ing a matrix of this size in computer memory requires 
at least 112Gb of RA M ( in single precision). One can 
follow the approach of [24| where we have computed the 
product K*K iteratively by storing only sufficiently small 
blocks of K in computer memory. Then the product 
K*K was Tikhonov-regularized by applying the opera¬ 
tion K*K —>■ K*K -I- \^I and the resultant system of 
equations was solved by the conjugate gradient descent. 
The bottleneck of this approach is the computation of 
K*K, which can still take very considerable time. 

On the other hand, solving the same linearized problem 
by the method described in Appendix B of [H is a much 
simpler task, and this approach yields, essentially, the 
same result. Indeed, in all cases we have considered, 
the two results were visually indistinguishable. This fact 
illustrates the proposition that considering the matrices 
A and B separately rather than combining them into one 
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large matrix K is computationally advantageous even in 
the linear regime. One can understand this improvement 
as a result of data reduction, that is, defining a transform 
of the data that is smaller in size but still contains all 
essential information. 

All linearized reconstructions shown below were ob¬ 
tained by regularizing and solving the equation W\v) = 
li.'exp) (sec. 7 of i)- The positive-definite matrix W of 
this equation was Tikhonov-regularized by the operation 
W W + \^ I. The regularization parameter A was ad¬ 
justed manually to obtain the best linear reconstruction 
in each case considered. The equation was then solved 
by the conjugate-gradient descent method. We note that 
the computational complexity of computing W (which 
is exactly of the same size as K*K, that is, Ny x Ny) 
is much smaller than the numerical complexity of com¬ 
puting K* K. Indeed, computing W requires only the 
singular vectors of the much smaller matrices A and B. 
Therefore, the main computational bottleneck of the lin¬ 
earized inversion is removed in this approach. 

We have applied different linearization methods, in¬ 
cluding first Born, first Rytov and mean-field approxi¬ 
mations as described in Appendix A of [i| . The choice of 
a linearization method only affects the way in which the 
data are computed and not the matrices AT or IT. Exten¬ 
sive simulations have revealed that reconstructions based 
on first Rytov or mean-field approximations do not pro¬ 
vide any noticeable advantages when compared to first 
Born. First Born and first Rytov reconstructions are 
compared in Fig. [2] below for illustrative purposes, but 
in all other cases only first Born-based reconstructions 
are shown. 


F. Error measures 

To quantify the convergence of the method, we use 
normalized root mean square errors (error of the so¬ 
lution) and r] 4 , (error of the equation). These quantities 
are defined as 

< = LlReconstructed) _ ^(True)l ^ ( 3 ^^) 

^ ^ • (31b) 

Note that t]^ is computed after Step 1 of the streamlined 
iteration cycle with the use of Computational Shortcut 
2 or after Step 2 for the algorithm without the use of 
Computational Shortcut 2 (as defined in Q). The error 
r],p is computed after Step 4 of this algorithm. After 
Step 5, the T-matrix is fully data-compatible and the 
error in question is zero up to the numerical precision of 
the computer. It should be kept in mind that can be 
computed only if the target is known a priori, which in 
practical applications is almost never the case. However, 


r] 4 , can be computed even if the target is not known a 
priori. 


IV. NUMERICAL RESULTS 
A. Small target 

Reconstructed images of the small target after 900 
DCTMC iterations are shown in Fig. [5] for various de¬ 
grees of contrast. Images obtained by linearized inver¬ 
sion are also shown for comparison. It can be seen that 
the linearized image reconstruction methods break down 
between xo = 0-01 and xo = 0-1. First Rytov approx¬ 
imation does not provide a visible advantage over first 
Born, and the same is true for the mean-field approxi¬ 
mation (data not shown). However, DCTMC yields re¬ 
constructions that are close to the correct result for all 
values of the contrast used. At xo = 1-75, the DCTMC 
reconstruction starts to visibly break down. The rea¬ 
son for this is incorrect assignment of three noninteract¬ 
ing voxels, as can be seen from the figure with sufficient 
magnification. The vicinity of these three voxels is also 
reconstructed incorrectly. So the breakdown in this case 
is due to the ad hoc algorithm for assigning the noninter¬ 
acting voxels. This problem is not fatal for DCTMC as 
the algorithm can be adjusted. Note that it is possible 
to learn whether some noninteracting voxels have been 
assigned incorrectly by looking at the convergence curve 
for the error 77 ,^ (see Fig. [3] below). Doing so does not 
require any a priori knowledge of the target. 

The errors rj^ and for the small target are shown in 
Fig. [3] as functions of the iteration number, i. We start 
the discussion of convergence with the error of recon¬ 
struction, rj^. The curves r]^(i) all look very similar and 
are almost independent of xo^ except for xo = 1-75. The 
latter case will be discussed separately. For xo < l-fb, 
the function ri^{i) has three distinct convergence regimes, 
as described below. 

First, there is a region of slow convergence. While 
DCTMC is in this regime, the reconstructed values Xn 
(or an) are all very small. Therefore, it can be con¬ 
cluded that the initial iterations simply solve the lin¬ 
earized inverse problem by the iterative algorithm de¬ 
scribed in [I| (Richardson first-order iteration). Natu¬ 
rally, convergence of this iterative process is expected 
to be slow. In principle, this initial slow convergence 
regime can be completely avoided by solving the lin¬ 
earized problem directly or by a fast iterative method 
such as conjugate-gradient descent, and then using the 
result as the initial guess for DCTMC. As is explained 
in Sec. IHI E[ linearized inversion can be obtained rela¬ 
tively fast with the use of data reduction that is inherent 
in treating the matrices A and B separately rather than 
combining them into one large matrix K. The computa¬ 
tional advantage obtained by this approach is illustrated 
in Sec. IV Al below. 

Second, there is a region of fast convergence. As- 
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FIG. 2. Linear (top and middle rows, marked FB and FR) and nonlinear (bottom row, marked NL) reconstructions of the 
small target for different levels of contrast xo. The quantity shown by the color scale in each plot is X"/xo, where Xn is the 
reconstructed susceptibility of the n-th voxel (real under the assumptions used) and xo is the amplitude of the shape function. 
For the linearized reconstructions, FB denotes first Born approximation and FR denotes first Rytov approximation. 


signment of noninteracting voxels occurs in this range 
of iteration indexes. If the curves are viewed with suf¬ 
ficient magnification, it can be seen that they are not 
smooth but contain ’’kinks”, that is, points where the 
slope changes abruptly. This change of slope occurs when 
one or more noninteracting voxels are determined cor¬ 
rectly. 

Third, there is the last region in which the rate of con¬ 
vergence is lower than in the second region but still much 
larger than in the first, slow convergence region. Clearly, 
the error %(f) decreases in the third region according to 
an exponential law. In this convergence regime, all non¬ 
interacting voxels have been determined correctly and 
the algorithm improves its estimate of the amplitudes of 
the remaining interacting voxels. The error continues 
to decrease to some very small values. This means that 
the reconstructions can be made very precise. However, 
the exponential convergence can not continue indefinitely 
because the error 77 ^^ can not decrease to an arbitrarily 
small value; it is bounded from below either by the ill- 
posedness of the inverse problem or by round-off errors. 

The above discussion is valid for the contrasts xo < 
1.75. For Xo = 1-75, the pattern is quite different. 
The reason is that three noninteracting voxels have been 
assigned incorrectly in this case. As a result, in the 


first, slow convergence region, the function %(i) actu¬ 
ally increases. These jumps take place when noninter¬ 
acting voxels are assigned incorrectly (apparently, too 
early). Then there still exists the fast convergence re¬ 
gion, wherein many noninteracting voxels are assigned 
correctly. Finally, the third, exponential convergence re¬ 
gion does not exist for xo = 1-75. This is so because the 
incorrectly assigned noninteracting voxels set a relatively 
large lower bound for the error r]^. 

We next discuss the error rj^. This error can be com¬ 
puted even if the target is not known a priori. First, 
we note that the dependence r],p{i) can also be classified 
into three different regimes (slow, fast and exponential), 
just as it was done for %(*). However, unlike in the 
case of 77x(Oi the curves r] 4 ,{i) depend noticeably on the 
contrast, xo- The dependence is, of course, weak for 
small Xo- Thus, the two curves for xo = 0.00175 and 
Xo = 0.0175 are barely distinguishable. However, the 
curves for xo = 0.0175, xo = 0.175 and xo = 0.875 can be 
easily distinguished. The difference is most pronounced 
in the exponential convergence region. The exponent ap¬ 
pears to be the same but the overall factor depends on 
Xo- This dependence of r],p on xo is a manifestation of 
the nonlinearity of the inverse problem. Indeed, it can 
be easily shown that, in the linear regime, xo 0 , 77 ^> is 
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independent of xo- We note that the dependence on xo 
can also be visible in the slow convergence region of the 
iteration indexes. This does not contradict the previously 
made observation that, in the slow convergence regime, 
DCTMC solves the linearized problem iteratively. The 
reason is that the definition of 77 $ involves the data ma¬ 
trix which, in general, is not proportional to xo- 

It was mentioned above that the overall factor of the 
function rj<p{i) in the exponential convergence region de¬ 
pends on Xo- Initially, this factor increases with xo- More 
generally, the curve 77 ^( 7 ) for xo = 0.875 goes higher than 
the curve for xo = 0.175 and the curve for xo = 0.175 
goes higher than the curve for xo = 0.0175, etc. How¬ 
ever, at larger values of xoj this tendency is reversed. 
Thus, the curve for xo = l-’75 starts lower than the curve 
for Xo = 0.00175 even before any noninteracting voxels 
have been assigned. This non-monotonous dependence 
oi ri<p on xo is a manifestation of the rather complicated 
nonlinearity of the inverse problem and is suggestive of 
a resonance phenomenon. That is, when we increase the 
value of xoj it passes close to a pole in the complex plane 
where the T-matrix has a singularity. We note that this 
non-monotonous dependence was observed for the large 
target as well (see below). 

The curves r]${i) are rather interesting in the strong 
nonlinearity case xo = 1.75. The error x <5 initially in¬ 
creases due to the incorrect assignment of the three non¬ 
interacting voxels. Then the error drops rapidly when 
many noninteracting voxels are assigned correctly. Then 
the error undergoes exponential decay with the same ex¬ 
ponent as for the smaller contrasts. Interestingly, the 
error 77 ^ in this range of iteration indexes is nearly con¬ 
stant while the error r],p steadily decreases. This is so 
because is dominated at this point by the three incor¬ 
rectly assigned noninteracting voxels while the algorithm 
still improves the accuracy of Xn in the remaining ma¬ 
jority of voxels. Finally, the exponential decay of 
crosses over to exponential growth. Why this happens is 
not entirely clear; one could expect r] 4 ,(i) to plateau. We 
can conjecture that this crossover to exponential growth 
occurs due to a complex interplay of the physical con¬ 
straint (that is still being applied at each iteration) and 
the use of an incomplete computational domain due to 
the incorrect exclusion of some of the voxels. We note in 
passing that the best reconstruction result for xo = l-’75 
would have been obtained if we stopped the iterations at 
i m 600. 

We conclude that the characteristic behavior of r] 4 ,{i) 
gives one an unambiguous indication that noninteracting 
voxels have been assigned incorrectly or (in more severe 
cases) that the obtained DCTMC reconstruction is not 
useful. 


B. Large target 

We now turn to the large target. In Fig. SI we show the 
reconstructions for five different values of the contrast. 



1 200 400 600 800 


FIG. 3. Convergence data for the small target. Errors and 
77 # are plotted vs the iteration number i. Different line types 
correspond to different contrast xo as follows: xo = 0.00175 
(1), xo = 0.0175 (2), xo = 0.175 (3), xo = 0.875 (4), and 
XO = 1.75 (5). 


from Xo = 0.002 to xo = 2. Of course, the computa¬ 
tional domain and the size of the inhomogeneities are in 
this case larger than in the small target and we can ex¬ 
pect the onset of nonlinearity to occur at smaller values 
of the contrast. Also, the linearized inverse problem is 
more ill-posed because there are voxels in the interior of 
the large target that are quite far from any source or de¬ 
tector and not effectively probed by incident evanescent 
waves. Indeed, the linearized reconstructions of the large 
target are not as good as the similar reconstructions of 
the small target. We emphasize that these are the best 
linearized reconstructions we were able to obtain by tun¬ 
ing the Tikhonov regularization parameter. Also, and 
as was the case for the small target, neither first Rytov 
nor the mean-field approximation provide a noticeable 
improvement of the linearized reconstructions (data not 
shown). 

We note however that the DCTMC reconstructions of 
the large target in the weak nonlinearity regime (e.g., 
Xo = 0.002) are considerably better than the linearized 
reconstructions. Moreover, it is not possible to improve 
the image quality of the linearized reconstruction by set¬ 
ting to zero the amplitudes of all voxels that are less in 
magnitude than, say, 1/40 of the maximum (recall that 
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40 is the smallest thresholding factor used by us in this 
paper to determine noninteracting voxels in DCTMC). 
In fact, a voxel with the reconstructed susceptibility 
Xn = Xo/40 is not visually distinguishable from zero in 
the color scheme used in this paper. Therefore, the differ¬ 
ence in quality between DCTMC and linearized methods 
is not a trivial consequence of image ” roughening”. The 
result may appear counter-intuitive. Indeed, in the limit 
Xo 0, the linearized inversion methods and DCTMC 
are solving exactly the same problem. However, the two 
approaches involve different regularization methods. In 
linearized inversions, Tikhonov regularization is used. In 
DCTMC, regularization is afforded by applying the phys¬ 
ical constraint at each iteration. Apparently, the latter 
approach is much better at reproducing sharp edges. 

As we move to the strong nonlinearity regime, the lin¬ 
earized reconstructions break down. At xo = 0 . 2 , the 
linearized reconstruction is not useful while DCTMC still 
provides a quantitatively accurate result. The contrast 
level Xo = 1 is borderline for DCTMC. It can be seen 
that the smaller, higher contrast inhomogeneity is still 
reconstructed correctly. However, the interior region of 
the larger inhomogeneity is not properly reconstructed 
(for the most part, underestimated). The boundaries of 
the larger inhomogeneity are nevertheless clearly visible. 
Note that we have encountered a similar situation in the 
case of nonlinear inversion by inverse Born series [^ . 
As in the case of the small target, the reason for the in¬ 
correct reconstruction of the larger inhomogeneity is the 
incorrect assignment of noninteracting voxels. However, 
the reconstruction is not yet entirely broken for the con¬ 
trast Xo = 1- At Xo = 2, too many noninteracting voxels 
have been assigned incorrectly and the resulting recon¬ 
struction is not useful. 

Convergence data for the large target are shown in 
Fig. [S] The behavior of the errors is qualitatively similar 
to what was observed in the case of the small target, al¬ 
though the overall rate of convergence is obviously lower. 
We note that the results for xo = 1 and Xo = 2 are qual¬ 
itatively different from those for smaller values of the 
contrast. This is explained by the incorrect assignment 
of noninteracting voxels at the higher levels of contrast. 

We note that the curve r]^{i) for the borderline case 
Xo = 1 has a well pronounced minimum at * = 590. 
We have observed a similar non-monotonous convergence 
for the smaller target as well. We have compared the 
reconstructions of the large target with Xo = 1 obtained 
at the “optimal” iteration index i = 590 and at i = 900 
9data not shown). As expected, the result is visibly but 
not dramatically better at i = 590. This result confirms 
our conjecture that monitoring the error is a useful 
approach to deciding when the iterations should stop. 


V. METHODS TO IMPROVE CONVERGENCE 

This section investigates several improvements to the 
“standard” DCTMC algorithm. It will be shown that. 


with a few simple modifications, we can substantially re¬ 
duce the number of iterations required to obtain accurate 
reconstructions. All simulations shown below concern the 
small target. 

A. Starting from the linear reconstrnction 

As was shown in [I|, DCTMC in the linear regime is 
equivalent to Richardson first-order iteration, which can 
have slow convergence. The regions of initial slow con¬ 
vergence are shown in Figs. [3] and [S] Two approaches can 
be used to cope with this problem. The first approach is 
to regularize the ISP as is explained in [l|. This, how¬ 
ever, can result in image degradation beyond what is war¬ 
ranted by the intrinsic ill-posedness of the problem. The 
second approach is to start iterations with the linearized 
solution to the ISP, V|in, as the initial guess. This is com¬ 
putationally efficient because we have a fast linearized 
solver at our disposal as is explained in Sec. IHIEI 

In Fig. [SJ we compare two initial guesses. The first 
initial guess is obtained by Vi = D[(/ -|- TexpF)“^T'exp]. 
The second initial guess is Vi = Idin- Both initial guesses 
were obtained for the small target and xo = 0.175. It can 
be seen that the first initial guess is very close to zero. 
It takes many Richardson iterations to transform it to 
the image resembling the one shown in the right panel. 
After this is accomplished, DCTMC starts to actually 
solve the nonlinear ISP. But the initial computational 
effort to arrive from the image in the left panel to the 
one in the right panel is unnecessary and can be avoided. 

To show that the number of necessary DCTMC iter¬ 
ations can be significantly reduced by using the initial 
guess Vi = Min, we have used the following algorithm: 

1: Obtain the linearized reconstruction. 

2: Using the result of Step 1 as the initial guess, run 
5 iterations of DCTMC normally. 

3: Then every 5 iterations check whether some of the 
values of Xn satisfy |xn| < Xmax/500, where Xmax = 
max„ |x„|. 

4: If a given voxel satisfies the above condition 3 
checks in a row, the corresponding Xn is set to zero, 
and the computational domain is reduced. 

5: The process is repeated with the following modifi¬ 
cations. After 15 iterations, checks are made every 
20 iterations and the relative threshold for deter¬ 
mining a non-interacting voxel is reduced to the 
factor of 100. After 200 iterations, checks are made 
every 10 iterations, and after 400 iterations, the rel¬ 
ative threshold is reduced to the factor of 60, and 
after 600 iterations the relative threshold is further 
reduced to the factor of 40. 

Compared to the original algorithm described above, 
the interval between the sparsity checks is reduced at 
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FIG. 4. Same as in Fig. O but for the large target and a somewhat different set of contrasts xo. Utilization of first Ry- 
tov approximation for linearized reconstruction does not provide any improvements over first Born approximation, and the 
corresponding results are not shown in this figure. 


the beginning to utilize the non-interacting voxels found 
by the linear reconstruction, while simultaneously rais¬ 
ing the relative threshold so as to not trust the linear re¬ 
construction too much. Empirical evidence suggests that 
linearized reconstructions tend to capture the boundaries 
of an objects more or less correctly but the interior and 
the regions close to the boundaries can be reconstructed 
with a large error 2^. Therefore, it can be expected 
that a linearized reconstruction would predict correctly 
the existing non-interacting voxels outside of the object 
and sufficiently far from its boundaries. However, nonin¬ 
teracting voxels in the interior or close to the boundaries 
can be predicted incorrectly. Therefore, we have modi¬ 
fied the first 15 iterations of the algorithm to avoid such 
incorrect assignments. After the initial 15 iterations, the 
parameters of the algorithm are reset to match the orig¬ 
inal method. 


The error plots illustrating the convergence of DCTMC 
when started from the two different initial guesses are 
shown in Fig. 0 A few immediate conclusions can be 
made from these plots. First, as we are starting from 
a much more reasonable guess, our initial error is much 
smaller. More importantly, as we can use sparsity checks 


accurately and relatively early in the iteration process, 
we have fast convergence almost from the start. Overall, 
for the contrast considered, we arrive at the final result 
in about 150 iterations less. This is significant but not 
yet dramatic. 

Also, the above improvement is not obtained for the 
strongest contrast we have considered, xo = 1.75. Here, 
when we start from the linearized reconstruction, we do 
not outperform the original method (data not shown). 
This can be understood qualitatively by noting that the 
linearized reconstruction at this contrast is very far from 
the true target. In fact, it is further away in the L 2 
norm from the true target than the zero initial guess. 
Moreover, the linearized reconstruction is very strongly 
affected by the choice of the Tikhonov regularization pa¬ 
rameter; an attempt to tune the latter to receive the 
best or most reasonable image quality can be counter¬ 
productive. More research is needed to understand how 
the best initial guess should be formed in the very strong 
nonlinearity regime. However, we will see that combin¬ 
ing all the improvements discussed in this section in one 
algorithm makes DCTMC convergence fast and reliable 
even at this high value of the contrast. 
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FIG. 5. Convergence data for the large target. The different 
line types correspond to different contrast xo as follows: xo = 
0.002 (1), xo = 0.02 (2), xo = 0.2 (3), xo = 1 (4) and xo = 2 
(5). 
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FIG. 6. Two different initial guesses (DCTMC iteration start¬ 
ing points) for the small target with xo = 0.175. The left 
panel is the initial guess for the T-matrix, Ti = Texp; the im¬ 
age is then computed by Vi = T>\{1 + rir)“^ri]. The right 
panel is Vi = Vun, where the linearized reconstruction Vin 
was obtained by using the first Born approximation. 


B. Using Reciprocity of Sources and Detectors 

For a vast majority of physical equations of interest, 
all measurements are reciprocal. In particular, this is 
true for the scalar wave equation considered here. This 
means that interchanging the source and detector yields 
the same measurement. For this reason, making physical 
measurements with interchanged sources and detectors 
does not provide any additional information about the 
target. Mathematically, reciprocity of measurements is 
equivalent to symmetry of the T-matrix T. 


FIG. 7. Gonvergence rate for the small target with xo = 
0.175. Different curves correspond to using the initial guess 
shown in Fig. [O] as labeled. 


In spite of what has been said above, there exists a 
subtle point in DCTMC that makes explicitly account¬ 
ing for the reciprocity useful. Namely, the experimen¬ 
tal T-matrix Texp is generally non-symmetric. There¬ 
fore, it does not account for the reciprocity of measure¬ 
ments. Consequently, including the data points with in¬ 
terchanged sources and detectors in the data matrix <I> 
will change Texp and make it closer to being symmetric 
and closer to correspondence to a diagonal V. This ac¬ 
count of reciprocity does not require any additional phys¬ 
ical measurements. All that is needed is to include in <P 
the data points that correspond to each source-detector 
pair being interchanged. In the traditional approaches to 
ISPs, such additional data points can be viewed as com¬ 
pletely redundant and useless. However, in the context 
of DCTMC, inclusion of these data points makes Texp 
more informative. In a sense, by including these redun¬ 
dant data points, we apply the additional condition that 
Texp should be as symmetric as possible in the real space 
representation. 

We now investigate the doubling of the data set that 
is obtained by interchanging sources and detectors. In 
fact, the consequences of this operation can be quite dra¬ 
matic. To illustrate the point, we show in Fig. [8] the con¬ 
vergence plots for the small target at the contrast level 
Xo = 0.0175. The initial guess Ti = Texp was used as the 
starting point. Two different data sets were used in the 
reconstructions. In the first case, sources are on one side 
of the sample and detectors are on the other side. In the 
second case, we supplement this data set by the “recip¬ 
rocal” data points obtained by interchanging the source 
and the detector in each source-detector pair. The size of 
the data set is doubled by this procedure. We emphasize 




































































13 



1 200 400 600 800 


FIG. 8. Convergence data for the small target with xo = 
0.0175 comparing the original DCTMC algorithm and the one 
using reciprocity of sources and detectors. 

again that the addition of the “reciprocal” data points is 
not truly redundant in DCTMC. Indeed, the fraction of 
known (non-zero) entries in the experimental T-matrix 
Texp (in the singular function representation) is 4% and 
17% in the two cases, respectively. That is, by using reci¬ 
procity, the dimension of the non-zero minor in Tgxp that 
is used in overwriting the iteratively updated T-matrix 
is roughly doubled. In other words, the experimental 
T-matrix is more informative in the second case. 

It can be seen that the convergence speed is much 
higher with the use of reciprocity. Moreover, the con¬ 
verged results are much more precise with both errors 
rjx and rj,p approaching the limit determined by the nu¬ 
merical precision of the computer. This is a dramatic 
improvement for a simple change that does not require 
any additional information or processing power. The case 
represented in the Figure is for xo = 0.0175, but all other 
levels of contrast that we have considered exhibit a sim¬ 
ilar behavior. 


C. Using weighted summation to the diagonal for 
the operator X>[ ] 

Next, we abandon the use of the Computational Short¬ 
cut 2 (defined in Ref. i) and use instead the force- 
diagonalization operator T>[-], viz, 

. (32) 

k 

Correspondingly, we will use the iteration scheme “Main 
iteration without the use of Computational Shortcut 
2,” [l|. This algorithm is more in line with the nonlocal 


framework of DCTMC. Although the computational time 
per one iteration will be increased by approximately the 
factor of 2, we will see that the improvements obtained 
by this approach are the most promising in terms of sig¬ 
nificantly improving the convergence rate of DCTMC. It 
is worth emphasizing that the results in this sub-section 
are obtained without the improvements related to start¬ 
ing from a better initial guesses and from doubling the 
data set by interchanging sources and detectors. We will 
also not apply the sparsity checks. Therefore, the results 
shown below can be attributed to DCTMC in its most 
basic form. 

As an example, we use a simple step-function for p{£). 
That is, p{£) = liti < R and p(^) = 0 otherwise. For the 
contrast xo = 0.0175, 150 DCTMC iterations were per¬ 
formed without any sparsity checks. The reconstructions 
obtained after these 150 iterations for varying values of R 
are shown in Fig. O The immediate conclusion is that the 
case R = 5 is much better than the original case R — 0. 
It is also true that R = lOh is too large, coupling vox¬ 
els that are too far apart, which creates a characteristic 
aliasing. 

The convergence data (for only and up to the it¬ 
eration index i = 150) are shown in Fig. [TO] It can be 
seen that a dramatic improvement of the convergence 
speed is achieved by using weighted summation to the 
diagonal with the appropriate radius of influence. The 
case i? = 5 is optimal, with R = 2.5 and R — 7.5 be¬ 
ing close in performance. Furthermore, it is interesting 
that for the case R = 10, the first couple of iterations 
are quite effective, before something goes wrong and the 
error curve turns sharply upwards. This example is not 
suppressing as the nonlocal interaction of the far away 
voxels causes an unwanted behavior. Note that the case 
R = 0 corresponds to the “standard” algorithm and the 
Computational Shortcut 2 was used for this case. But 
despite the fact that the R = 0 case completed in about 
half the time, looking at the error plot in Fig. [TUI shows 
that there is no confusion on the preference between 75 
iterations for i? = 5 and 150 iterations of i? = 0. 

D. Putting it all together 

We now combining all three improvements described 
above in one algorithm. We will demonstrate that for all 
values of the contrast used for the small sample, we can 
achieve acceptable either perfect or at least acceptable 
results in 75 iterations, compared to the 900 iterations 
previously performed. The algorithm used in this section 
proceeds as follows: 

1: Run the linear reconstruction. 

2: Using the result of step 1 as the initial guess, run 
5 iterations without sparsity checks. 

3: Then every 5 iterations check whether some sus¬ 
ceptibilities Xn satisfy Ixnl < Xmax/500, where 
Xmax = max„ \x„\. 
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R = 0 R — 3h R — 5h R — 7h R = lOh 


FIG. 9. DCTMC reconstruction of the small target after 150 iterations with varying degrees of R, where the row-summing is 
done over voxels separated by no more than R. 



FIG. 10. Convergence data for the reconstruction shown 
in Fig. m 


4: If a given voxel satisfies the above condition 3 
checks in a row, the corresponding Xn is set to zero, 
and the computational domain is reduced. 

5: The process is repeated with the following modifi¬ 
cations. After 20 iterations, the relative threshold 
for determining a non-interacting voxel is reduced 
to the factor of 100. After 40 iterations, the relative 
threshold is further reduced to the factor of 60. 

Throughout the algorithm, reciprocity of sources and 
detectors was used, as well as weighted summation to 
the diagonal for the operator 'D[-]. For the weight func¬ 
tion, we used the inverse distanced, namely, p{£ik) = 

+ (1 ~ Computational Shortcut 2 was not 

used. The reconstruction after 75 iterations are shown 
in Fig. [TT] for all levels of contrast used previously. All 
reconstructions are nearly perfect at least for the largest 
value of the contrast, xo = 1-75. In the latter case, the 
reconstruction is not perfect but acceptable. The impre¬ 
cision is again caused by three voxels that were incor¬ 
rectly assigned as noninteracting. This happened even 
though we have increased the minimum relative thresh¬ 
old for determining noninteracting voxels from 40 to 60. 

The convergence data for the algorithm employed in 
this subsection are shown in Fig. [12] It can be seen 
that the most dramatic improvement was obtained for 
the three lowest levels of contrast, as starting from the 
linearized inversion is of great assistance in these cases. 


In the intermediate case xo = 0.875, the reconstruction 
after 75 iterations is not as precise as the one produced by 
the original method after 900 iterations. This is clearly so 
because one voxel was erroneously assigned as noninter¬ 
acting by the improved method. One incorrectly assigned 
noninteracting voxel does not degrade the reconstructed 
image seriously but bounds from below the error 77 ^ that 
can be reached. Clearly, more research on the applica¬ 
tion of sparsity checks in the strong nonlinearity regime 
is needed, especially when linearized inversion is used as 
the initial guess.. 

VI. COMPARISON WITH GAUSS-NEWTON 
METHOD 

For comparison purposes, we have also implemented 
the regularized Gauss-Newton method. We used the 
standard implementation of Levenberg-Marquardt with 
multiplicative damping. This method was too slow to 
be applied to the large target but we could apply it to 
the small target, even though one iteration of the Gauss- 
Newton method took much longer than one iteration of 
DCTMC. More importantly, the Gauss-Newton method 
did not converge to a reasonable result for the small tar¬ 
get even at the lowest level of the contrast considered 
above. It is worth stressing that numerous modifications 
of the regularization and damping techniques were tried 
without producing any noticeably better results. There¬ 
fore, the Gauss-Newton reconstruction shown below is 
representative. 

A typical Gauss-Newton reconstruction of the small 
target with ap = 0.00175 is shown in Fig.[T31 The corre¬ 
sponding error plots are shown in Fig. 1141 It is obvious 
from the data of Fig.[T4|that the Gauss-Newton iterations 
are converging to a false solution (a local minimum of the 
cost function). Specifically, the error of the reconstruc¬ 
tion rj^ is increasing with each iteration while the error 
of the equation r]$ is decreasing and reaching some very 
small values that are consistent with the limit imposed 
by the finite numerical precision of the computer. 

To confirm that the Gauss-Newton method was pro¬ 
grammed correctly, we have reduced the size of the tar¬ 
get even further to a point where the method works well, 
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FIG. 11. Reconstructions of the small target by the algorithm combining all three improvements discussed in this section after 
75 iterations for the five levels of contrast previously considered. 
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FIG. 12. Convergence data of for the algorithm combining 
all improvements discussed in this section. The horizontal 
lines show the respective error attained after 900 iterations 
of DCTMC without any of the improvements. Different line 
types correspond to different contrast xo as follows: xo = 
0.00175 (1), xo = 0.0175 (2), xo = 0.175 (3), xo = 0.875 (4), 
and Xo = 1.75 (5). 


at least, for some levels of contrast. The shape of this 
even smaller (tiny) target is shown in Fig. [T] and recon¬ 
structions are shown in Fig. [15] For the tiny target, the 
Gauss-Newton reconstructions work well except at the 
high level of contrast (xo = 3), when the method fails to 
localize the inhomogeneity correctly. 

Thus, we have produced numerical evidence that the 
DCTMC algorithm is able to solve nonlinear ISPs where 
other mainstream nonlinear inversion techniques fail. 
The results reported here further suggest that the size 
of the target is a critical factor. Specifically, the Gauss- 
Newton method works reasonably well for the tiny tar¬ 
get but fails dramatically for the small target. A pos- 


FIG. 13. A typical Gauss-Newton (with Levenburg- 
Marquardt damping technique) reconstruction of the small 
target with oq = 0.00175. 



FIG. 14. Convergence data for the Gauss-Newton (with 
Levenberg-Marquardt damping technique) reconstruction 
shown in Fig. 1131 Note that for i > 5, the error rj<p{i) contin¬ 
ues to decrease, as is guaranteed by the algorithm, although 
this decrease is not visible due to the logarithmic vertical scale 
used in the figure. 

















































































































































































































































16 



FIG. 15. Gauss-Newton (with Levenburg-Marquardt damp¬ 
ing technique) reconstructions of the tiny target for different 
levels of contrast as labeled. 


sible mathematical interpretation of this empirical ob¬ 
servation is that the cost function used by the Gauss- 
Newton method develops a spurious local minimum for 
the small target. This minimum is connected by some 
path in the multi-dimensional space of solutions along 
which the cost function increases to the linearized solu¬ 
tion. It follows therefore that the linearized solution is 
not a true minimum of the cost function. This state¬ 
ment may seem counterintuitive since we know for sure 
that the linearized solution minimizes the linearized cost 
function. However, the subtle but important point is 
that we do not really know whether the linearized solu¬ 
tion minimizes the nonlinear cost function. It is also a 
realistic possibility that the linearized solution is in fact 
a minimum of the nonlinear cost function but this mini¬ 
mum is so shallow that, numerically, the iterations always 
overshoot and eventually descend towards the false solu¬ 
tion. Clearly, more research in this direction is needed 
to fully understand the advantages and disadvantages of 
DCTMC compared to the standard methods. 


VII. DISCUSSION 


We have provided an initial numerical investigation 
of the data-compatible T-matrix completion (DCTMC) 
method for solving nonlinear inverse scattering problems 
(ISPs). DCTMC was implemented for the scalar wave 
equation of the form ([2]) which is relevant, in particular, 
to ultrasound imaging and seismic tomography. In many 
of the cases we have considered, DCTMC provided quan¬ 
titatively accurate reconstructions whilst any linearized 
inversion (including those based on first Born, first Ry- 
tov or mean-field approximations) or the standard op¬ 
timization methods (Gauss-Newton) have failed. This 
provides us with a proof-of-principle demonstration of 
the method’s utility. Indeed, DCTMC does not utilize a 
cost function and does not suffer from the local minima 
and the associated false solution. Perhaps, this is the 
most important conceptual difference that distinguishes 
DCTMC from other nonlinear inversion methods. The 
method is especially useful for solving ISPs with large 
data sets. Even though we did not use supercomput¬ 


ers or any massively parallel computational platforms, 
we have succeeded in solving a strongly nonlinear inverse 
problem with more than 2 • 10® data points. Still, many 
questions regarding DCTMC improvement, optimization 
and applicability remain to be investigated either theo¬ 
retically or numerically. Some of these topics for future 
research are listed below. 

First, the obvious task is to investigate DCTMC with 
adjustable sparsity checks. These checks are very use¬ 
ful for improving the speed and convergence of DCTMC. 
However, incorrect assignment of noninteracting voxels 
has resulted in the method breakup in the strong nonlin¬ 
earity regime. We emphasize that all cases of DCTMC 
not providing a useful reconstruction were due to these 
incorrect noninteracting voxel assignments. If the spar¬ 
sity checks are suppressed, the method can be slowly con¬ 
verging but we have not encountered and divergences or 
convergence to false solutions. 

Second, we did not investigate in much detail various 
approaches to regularization. It is obvious that, as the 
target size increases, the ISP becomes more ill-posed. We 
can attempt a combination of Tikhonov regularization of 
the type that was discussed in [I[, physical constraints 
and, potentially, other methods. However, complicated 
interaction between different types of regularization re¬ 
quire a careful and systematic investigation. 

Third, there are questions of computational efficiency 
and inclusion of larger targets. A combination of meth¬ 
ods described in Sec. 13 is expected to increase the rate 
of convergence or reduce the number of necessary iter¬ 
ations. Improving the computational bottlenecks is an¬ 
other important goal for future research. We note that 
some approximate approaches to speeding up this aera¬ 
tions have been proposed by Jakobsen and Ursin @ and 
these approaches can be utilized by DCTMC. 

Fourth, weight functions p{t) used in the definition of 
the force-diagonalization operator 'D[-] require a separate 
investigation 

We finally note that the major challenge in understand¬ 
ing and optimizing the DCTMC is the complicated inter¬ 
play of the various tunable parameters and algorityhms 
thart are inherent to the method. Nevertheless, the au¬ 
thors are optimistic that further progress can be made by 
a combination of analytical and numerical investigations. 
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